Sapling Linear Regression
Kathleen Stutzman
2024-01-10
Libraries:
# libraries -------------------
library(dplyr)
library(tidyverse)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(TMB)
library(emmeans)
library(DHARMa)
library(lmtest)
library(ggeffects)
library(performance)Set working directory
knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')Source sapling data:
source("Scripts/SA_Import.R")Pivot wider to create dataframe where each row is for one plot and has total details for each species group
sa_merge2 <- sapling_merge %>%
pivot_wider(names_from = Species_Groups, values_from = Total_Tally)Import time since data and add it to the PIRI dataset
time_since <- read_csv("CleanData/Treat_Year_Data.csv")
sa_merge3 <- merge(sa_merge2, time_since, by = 'Site')
#log transform time from treatment data
sa_merge3$l.TFT <- log(sa_merge3$Time_from_Treat)Run the ‘Add_BA’ script and merge with dataset:
source("Scripts/Add_BA.R")
# merge with sapling dataset -------------------
sa_merge4 <- merge(sa_merge3, prism_BA, by = 'Plot_No')Run ‘Ground_Data.R’ script and add it to sapling dataset:
source("Scripts/Ground_Data.R")
# merge with sapling dataset -------------------
sa_merge5 <- merge(sa_merge4, ground3, by = 'Plot_No')Source and add veg data:
source("Scripts/Veg_Data.R")
# merge with sapling dataset -------------------
sa.all <- merge(sa_merge5, veg3, by = "Plot_No")
rm(sa_merge3,
sa_merge4,
sa_merge5)Sapling count data is taken in 25m2 plots; basal area is measured in hectares; veg and soil data is taken a 1m2 plots.
Data will be standardized for comparisons at multiple scales
Create a dataframe keeping sapling at 25m2 observation levels & log transforming them
sa.all3 <- sa.all
sa.all3$piri.ba_s <- scale(sa.all3$PIRI.BA_HA)
sa.all3$avgld_s <- scale(sa.all3$avgLD)
sa.all3$vegt_s <- scale(sa.all3$Veg_Total)
sa.all3$min_s <- scale(sa.all3$Mineral_Soil)
sa.all3$so_s <- scale(sa.all3$Shrub_Oak)
sa.all3$other_s <- scale(sa.all3$Other)
sa.all3$ba.ha_s <- scale(sa.all3$BA_HA)
sa.all3$piri_s <- scale(sa.all3$PIRI)
sa.all3 <- sa.all3 %>%
select(Treat_Type, Region, Site, Plot_No, PIRI, piri_s, Shrub_Oak, so_s, Other, other_s, Time_from_Treat, l.TFT, BA_HA, ba.ha_s, PIRI.BA_HA, piri.ba_s, Mineral_Soil, min_s, Litter_Duff, avgLD, avgld_s, Veg_Total, vegt_s) %>%
arrange(Treat_Type)Looking at pairwise with sa.all3 dataset (sapling counts at 25m2 observation level)
#not log transformed (25m2)
sa3.num <- sa.all3 %>%
select(PIRI, Shrub_Oak, Other, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)
ggpairs(sa3.num)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggpairs(sa3.num, aes(color = Treat_Type))## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
rm(sa3.num)No strong correlations
major personal breakthrough for dummies. For zero inflation, lets see if any TT or region is lacking alltogether in PIRI or SO saplings
sa.all4 <- sa.all3 %>%
group_by(Region)
tapply(sa.all4$PIRI, sa.all4$Region, summary)## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.02459 0.00000 2.00000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2013 0.0000 9.0000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.07273 0.00000 1.00000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.09574 0.00000 3.00000# no PIRI in LI,
tapply(sa.all4$PIRI, sa.all4$Treat_Type, summary) ## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2564 0.0000 9.0000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.04902 0.00000 2.00000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.05357 0.00000 2.00000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.05882 0.00000 2.00000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01149 0.00000 1.00000Now that I know that there are no PIRI in LI, I want to start model elimination again using ziformula = ~ Region
I have previously done model elimination with zero inflated poisson distro & with model w.o TT. Neither are better than zi nbinom2 models. zi ~1 has distribution issues, ~Region does not
Model has much better fit and zero inflation issues are resolved when treatment type variable is included. Time to do elimination again but with TT; all models fail with nbinom1 distro, but works with nbinom2 (these notes are from when I was using a zi poisson distro, using nbinom TT no longer seems to be important to model fit, but we will keep for the overarching story)
Checking fixed effects for collinearity and looking at general summary stats to make sure they aren’t 0/1 or same/function of other category
summary(sa.all3)## Treat_Type Region Site Plot_No
## Length:498 Length:498 Length:498 Min. : 70.0
## Class :character Class :character Class :character 1st Qu.: 359.5
## Mode :character Mode :character Mode :character Median : 614.0
## Mean : 613.5
## 3rd Qu.: 887.5
## Max. :1143.0
## PIRI piri_s.V1 Shrub_Oak so_s.V1
## Min. :0.00000 Min. :-0.157452 Min. :0.0000 Min. :-0.273816
## 1st Qu.:0.00000 1st Qu.:-0.157452 1st Qu.:0.0000 1st Qu.:-0.273816
## Median :0.00000 Median :-0.157452 Median :0.0000 Median :-0.273816
## Mean :0.09438 Mean : 0.000000 Mean :0.2691 Mean : 0.000000
## 3rd Qu.:0.00000 3rd Qu.:-0.157452 3rd Qu.:0.0000 3rd Qu.:-0.273816
## Max. :9.00000 Max. :14.857482 Max. :7.0000 Max. : 6.849481
## Other other_s.V1 Time_from_Treat l.TFT
## Min. :0.0000 Min. :-0.345712 Min. : 1.0 Min. :0.000
## 1st Qu.:0.0000 1st Qu.:-0.345712 1st Qu.: 3.0 1st Qu.:1.099
## Median :0.0000 Median :-0.345712 Median : 5.0 Median :1.609
## Mean :0.3514 Mean : 0.000000 Mean :10.9 Mean :1.804
## 3rd Qu.:0.0000 3rd Qu.:-0.345712 3rd Qu.: 8.0 3rd Qu.:2.079
## Max. :9.0000 Max. : 8.508476 Max. :33.0 Max. :3.497
## BA_HA ba.ha_s.V1 PIRI.BA_HA piri.ba_s.V1
## Min. : 0.00 Min. :-1.4237625 Min. : 0.0 Min. :-1.285691
## 1st Qu.: 7.00 1st Qu.:-0.8185636 1st Qu.: 5.0 1st Qu.:-0.775675
## Median :16.00 Median :-0.0404508 Median :11.0 Median :-0.163656
## Mean :16.47 Mean : 0.0000000 Mean :12.6 Mean : 0.000000
## 3rd Qu.:23.00 3rd Qu.: 0.5647481 3rd Qu.:18.0 3rd Qu.: 0.550367
## Max. :51.00 Max. : 2.9855436 Max. :44.0 Max. : 3.202451
## Mineral_Soil min_s.V1 Litter_Duff avgLD
## Min. : 0.0000 Min. :-0.197640 Min. : 0.00 Min. : 0.000
## 1st Qu.: 0.0000 1st Qu.:-0.197640 1st Qu.:90.50 1st Qu.: 2.250
## Median : 0.0000 Median :-0.197640 Median :90.50 Median : 3.875
## Mean : 0.7199 Mean : 0.000000 Mean :88.78 Mean : 4.179
## 3rd Qu.: 0.0000 3rd Qu.:-0.197640 3rd Qu.:90.50 3rd Qu.: 5.500
## Max. :50.5000 Max. :13.666913 Max. :90.50 Max. :13.250
## avgld_s.V1 Veg_Total vegt_s.V1
## Min. :-1.652607 Min. : 4.00 Min. :-1.722482
## 1st Qu.:-0.762793 1st Qu.: 37.50 1st Qu.:-0.738485
## Median :-0.120151 Median : 58.00 Median :-0.136337
## Mean : 0.000000 Mean : 62.64 Mean : 0.000000
## 3rd Qu.: 0.522492 3rd Qu.: 83.38 3rd Qu.: 0.609004
## Max. : 3.587405 Max. :202.00 Max. : 4.093383library(caret)## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:purrr':
##
## liftfindLinearCombos(sa.all3[,4:23]) #shows no linear combos. Don't get how the rank deficiency is happening## $linearCombos
## $linearCombos[[1]]
## [1] 5 2 3 4
##
## $linearCombos[[2]]
## [1] 7 2 3 6
##
## $linearCombos[[3]]
## [1] 11 2 3 10
##
## $linearCombos[[4]]
## [1] 13 2 3 12
##
## $linearCombos[[5]]
## [1] 15 2 3 14
##
## $linearCombos[[6]]
## [1] 18 2 3 17
##
## $linearCombos[[7]]
## [1] 20 2 3 19
##
##
## $remove
## [1] 5 7 11 13 15 18 20tapply(sa.all3$Other, sa.all3$Region, summary) #all present## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.3443 0.0000 6.0000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.411 0.000 6.000
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2143 0.0000 6.0000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.6727 0.0000 9.0000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.3511 0.0000 4.0000tapply(sa.all3$Shrub_Oak, sa.all3$Region, summary) #none in LI## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01639 0.00000 1.00000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.000000 0.000000 0.006494 0.000000 1.000000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 1.964 3.000 7.000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2447 0.0000 6.0000tapply(sa.all3$BA_HA, sa.all3$Region, summary) #all present## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 11.00 14.98 23.00 51.00
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 11.00 21.00 19.81 28.00 39.00
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 9.00 14.00 15.32 21.00 44.00
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 13.50 23.00 22.13 31.00 44.00
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 11.00 14.37 21.00 48.00tapply(sa.all3$PIRI.BA_HA, sa.all3$Region, summary) #all present## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 2 9 10 16 41
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 11.00 16.00 17.41 23.00 37.00
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 9.00 11.44 16.00 34.00
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 8.00 18.00 19.18 30.00 44.00
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 9.00 10.32 16.00 32.00tapply(sa.all3$Veg_Total, sa.all3$Region, summary) #all present## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.00 51.50 71.00 70.96 91.88 166.50
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.00 23.00 39.00 44.52 61.00 128.00
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.50 33.75 52.75 56.93 74.38 161.50
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 14.50 41.25 64.00 62.97 79.75 114.00
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.00 46.00 62.25 75.09 96.62 202.00tapply(sa.all3$Mineral_Soil, sa.all3$Region, summary) #all present## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 1.586 0.500 30.500
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.589 0.500 15.500
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.6623 0.0000 50.5000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01818 0.00000 0.50000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2021 0.0000 8.0000tapply(sa.all3$avgLD, sa.all3$Region, summary) #all present## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.000 1.750 1.977 2.750 7.250
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.75 4.00 5.25 5.61 6.75 12.50
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.250 3.562 5.000 5.347 6.750 13.250
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.750 2.500 3.500 3.686 4.500 9.750
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.750 2.812 4.000 4.298 5.438 12.250tapply(sa.all3$Other, sa.all3$Treat_Type, summary) #all present## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.8889 1.0000 6.0000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.4118 0.0000 9.0000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.1429 0.0000 2.0000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.06618 0.00000 2.00000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.1379 0.0000 2.0000tapply(sa.all3$Shrub_Oak, sa.all3$Treat_Type, summary) #none in mowrx or springrx## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.7094 0.0000 7.0000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.4902 0.0000 6.0000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01786 0.00000 1.00000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0tapply(sa.all3$BA_HA, sa.all3$Treat_Type, summary) #all present## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 18.00 28.00 26.78 37.00 51.00
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 9.00 16.00 17.44 23.00 44.00
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 11.00 12.11 18.00 34.00
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.000 7.000 9.059 14.000 37.000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 10.00 14.00 15.85 21.00 41.00tapply(sa.all3$PIRI.BA_HA, sa.all3$Treat_Type, summary) #all present## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 9.00 18.00 17.05 25.00 41.00
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 11.00 14.45 20.25 44.00
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 5.000 9.000 9.804 14.500 30.000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.000 7.000 7.816 11.000 37.000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 9.00 14.00 13.75 18.00 34.00tapply(sa.all3$Veg_Total, sa.all3$Treat_Type, summary) #all present## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.00 26.50 46.00 46.85 61.00 115.50
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 14.5 46.0 68.0 72.7 91.0 202.0
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 15.00 36.12 55.00 55.21 68.62 132.00
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.00 37.00 66.75 70.47 94.38 196.50
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 16.00 49.50 61.50 64.65 81.00 131.00tapply(sa.all3$Mineral_Soil, sa.all3$Treat_Type, summary) #all present## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.000000 0.000000 0.008547 0.000000 0.500000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.201 0.000 8.000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 1.652 0.000 50.500
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 1.199 0.500 30.500
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.9368 0.5000 15.5000tapply(sa.all3$avgLD, sa.all3$Treat_Type, summary) #all present## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.500 3.250 4.750 5.024 6.750 12.500
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.750 3.250 4.375 4.529 5.750 9.750
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.250 3.688 4.775 5.291 6.750 11.000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.500 2.750 2.816 3.812 13.250
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.500 1.500 3.500 4.046 5.250 11.500PIRI Model
Best Model AIC: 245.8
Running model elimination used standardized effects and zi by Region
sa.m1s <- glmmTMB(PIRI ~ Treat_Type + piri.ba_s + avgld_s + vegt_s + min_s + so_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; false convergence (8). See vignette('troubleshooting'),
## help('diagnose')#won't converge
#test piri ba
sa.m2s <- glmmTMB(PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m2s) #248## [1] 247.957# test ba/ha
sa.m3s <- glmmTMB(PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m3s) #251.7## [1] 251.6733lrtest(sa.m2s, sa.m3s) #keep ba## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s +
## ba.ha_s + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 19 -104.98
## 2 18 -107.84 -1 5.7163 0.01681 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# test SO
sa.m4s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + other_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m4s) #248## [1] 248.2719lrtest(sa.m4s, sa.m2s) #drop## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + other_s +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s +
## ba.ha_s + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 18 -106.14
## 2 19 -104.98 1 2.3148 0.1281# test other
sa.m5s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m5s) #251.1## [1] 251.0924lrtest(sa.m4s, sa.m5s) #keep other## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + other_s +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 18 -106.14
## 2 17 -108.55 -1 4.8206 0.02812 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# test mineral
sa.m6s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + other_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m6s) #247## [1] 247.0116lrtest(sa.m4s, sa.m6s) #drop## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + other_s +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + other_s + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 18 -106.14
## 2 17 -106.51 -1 0.7397 0.3897# test veg
sa.m7s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + avgld_s + other_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m7s) #245.8## [1] 245.8436# test ld
sa.m8s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + other_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m8s) #244.8 ## [1] 244.7564sa.m8s_sr <- simulateResiduals(sa.m8s, n = 1000, plot = TRUE) #passes
testResiduals(sa.m8s_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.029239, p-value = 0.7881
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.4242, p-value = 0.826
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00072289156626506 )
## 0## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.029239, p-value = 0.7881
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.4242, p-value = 0.826
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00072289156626506 )
## 0# passes
emmeans(sa.m8s, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response') # the model is rank deficient## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.0272 NaN Inf NaN NaN
## FallRx 0.0235 NaN Inf NaN NaN
## Harvest 0.1676 NaN Inf NaN NaN
## MowRx 0.1509 NaN Inf NaN NaN
## SpringRx 0.0193 NaN Inf NaN NaN
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 1.157 NaN Inf 1 NaN NaN
## Control / Harvest 0.162 NaN Inf 1 NaN NaN
## Control / MowRx 0.180 NaN Inf 1 NaN NaN
## Control / SpringRx 1.409 NaN Inf 1 NaN NaN
## FallRx / Harvest 0.140 NaN Inf 1 NaN NaN
## FallRx / MowRx 0.156 NaN Inf 1 NaN NaN
## FallRx / SpringRx 1.217 NaN Inf 1 NaN NaN
## Harvest / MowRx 1.111 NaN Inf 1 NaN NaN
## Harvest / SpringRx 8.678 NaN Inf 1 NaN NaN
## MowRx / SpringRx 7.812 NaN Inf 1 NaN NaN
##
## P value adjustment: tukey method for comparing a family of 2 estimates
## Tests are performed on the log scalesa.m9s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + other_s + avgld_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m9s) #245.9 with veg; 245.8 with ld## [1] 245.8436summary(sa.m9s)## Family: nbinom2 ( log )
## Formula:
## PIRI ~ Treat_Type + ba.ha_s + other_s + avgld_s + offset(l.TFT) +
## (1 | Site/Plot_No)
## Zero inflation: ~Region
## Data: sa.all3
##
## AIC BIC logLik deviance df.resid
## 245.8 313.2 -106.9 213.8 482
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 0.00000003136 0.0001771
## Site (Intercept) 1.00553141354 1.0027619
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Dispersion parameter for nbinom2 family (): 0.332
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.4223 0.8811 -6.154 0.000000000755 ***
## Treat_TypeFallRx -0.2051 0.9401 -0.218 0.8273
## Treat_TypeHarvest 1.8988 1.3009 1.460 0.1444
## Treat_TypeMowRx 1.6541 1.0395 1.591 0.1116
## Treat_TypeSpringRx -0.2551 1.4500 -0.176 0.8603
## ba.ha_s 0.8721 0.3623 2.407 0.0161 *
## other_s -1.1511 0.6811 -1.690 0.0910 .
## avgld_s -0.2776 0.2909 -0.954 0.3398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.521 1.019 1.492 0.136
## RegionLI 22.406 53458.279 0.000 1.000
## RegionMA -20.529 22550.145 -0.001 0.999
## RegionME -2.795 3.323 -0.841 0.400
## RegionNH -1.687 1.457 -1.158 0.247Best model (lowest AIC without rank deficiency) appears to be model sa.m6 with TT, piri BA, agvLD, Veg total. Test model fit:
Model without avgld is rank deficient
sa.m9s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + other_s + avgld_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m9s) #245.9 with veg; 245.8 with ld## [1] 245.8436summary(sa.m9s)## Family: nbinom2 ( log )
## Formula:
## PIRI ~ Treat_Type + ba.ha_s + other_s + avgld_s + offset(l.TFT) +
## (1 | Site/Plot_No)
## Zero inflation: ~Region
## Data: sa.all3
##
## AIC BIC logLik deviance df.resid
## 245.8 313.2 -106.9 213.8 482
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 0.00000003136 0.0001771
## Site (Intercept) 1.00553141354 1.0027619
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Dispersion parameter for nbinom2 family (): 0.332
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.4223 0.8811 -6.154 0.000000000755 ***
## Treat_TypeFallRx -0.2051 0.9401 -0.218 0.8273
## Treat_TypeHarvest 1.8988 1.3009 1.460 0.1444
## Treat_TypeMowRx 1.6541 1.0395 1.591 0.1116
## Treat_TypeSpringRx -0.2551 1.4500 -0.176 0.8603
## ba.ha_s 0.8721 0.3623 2.407 0.0161 *
## other_s -1.1511 0.6811 -1.690 0.0910 .
## avgld_s -0.2776 0.2909 -0.954 0.3398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.521 1.019 1.492 0.136
## RegionLI 22.406 53458.279 0.000 1.000
## RegionMA -20.529 22550.145 -0.001 0.999
## RegionME -2.795 3.323 -0.841 0.400
## RegionNH -1.687 1.457 -1.158 0.247sa.m9s_sr <- simulateResiduals(sa.m9s, n = 1000, plot = TRUE) #passes
testResiduals(sa.m9s_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.030633, p-value = 0.7383
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.30212, p-value = 0.862
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000883534136546185 )
## 0## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.030633, p-value = 0.7383
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.30212, p-value = 0.862
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000883534136546185 )
## 0# passesPairwise
emmeans(sa.m9s, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.0268 0.0236 Inf 0.00477 0.151
## FallRx 0.0219 0.0206 Inf 0.00345 0.139
## Harvest 0.1792 0.1821 Inf 0.02444 1.314
## MowRx 0.1403 0.1040 Inf 0.03279 0.600
## SpringRx 0.0208 0.0280 Inf 0.00148 0.292
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 1.228 1.154 Inf 1 0.218 0.9995
## Control / Harvest 0.150 0.195 Inf 1 -1.460 0.5889
## Control / MowRx 0.191 0.199 Inf 1 -1.591 0.5030
## Control / SpringRx 1.291 1.871 Inf 1 0.176 0.9998
## FallRx / Harvest 0.122 0.166 Inf 1 -1.542 0.5352
## FallRx / MowRx 0.156 0.170 Inf 1 -1.704 0.4314
## FallRx / SpringRx 1.051 1.577 Inf 1 0.033 1.0000
## Harvest / MowRx 1.277 1.562 Inf 1 0.200 0.9996
## Harvest / SpringRx 8.619 14.264 Inf 1 1.302 0.6903
## MowRx / SpringRx 6.748 10.028 Inf 1 1.285 0.7008
##
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scaleNone significantly different from eachother
Plot model
sa.piri1 <- ggpredict(sa.m9s, terms = c("avgld_s", "Treat_Type"))
#control and fallrx are almost completely identical - if I wanted to nudge/jitter
# sa.piri2 <- sa.piri1 %>%
# filter(group == "FallRx")
#
# sa.piri2$predicted <- sa.piri2$predicted+0.01
#
# sa.piri3 <- sa.piri1 %>%
# filter(group != "FallRx")
#
# sa.piri3 <- rbind(sa.piri3, sa.piri2)
ggplot(sa.piri1, aes(x, predicted, color = group))+
geom_line(linewidth = 1, show.legend = FALSE) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(axis.text = element_text(size = 18), axis.title = element_text(size = 24))+
labs(x = "Average leaf litter depth \n (standardized)",
y = NULL)
ggsave(filename = "SA_PIRI_avgld.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)
sa.piri4 <- ggpredict(sa.m9s, terms = c("ba.ha_s", "Treat_Type"))
#again control and fall rx overlap - if I wanted to nudge/jitter
# sa.piri5 <- sa.piri4 %>%
# filter(group == "FallRx")
#
# sa.piri5$predicted <- sa.piri5$predicted+0.01
#
# sa.piri6 <- sa.piri4 %>%
# filter(group != "FallRx")
#
# sa.piri6 <- rbind(sa.piri6, sa.piri5)
ggplot(sa.piri4, aes(x, predicted, color = group))+
geom_line(linewidth = 1, show.legend = FALSE) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(axis.text = element_text(size = 18), axis.title = element_text(size = 24))+
labs(x = "Average basal area per hectare \n (standardized)",
y = "Pitch pine stems \n (adjusted for time)")
ggsave(filename = "SA_PIRI_BA.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)
sa.piri5 <- ggpredict(sa.m9s, terms = c("other_s", "Treat_Type"))
ggplot(sa.piri5, aes(x, predicted, color = group))+
geom_line(linewidth = 1, show.legend = FALSE) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(axis.text = element_text(size = 18), axis.title = element_text(size = 24))+
labs(x = "Counts of other saplings \n (standardized)",
y = NULL)
ggsave(filename = "SA_PIRI_other.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)Shrub Oak Models
Best Model AIC: 323.8
Start from the beginning knowing that both Region and Treat Type have 0 areas
tapply(sa.all4$Shrub_Oak, sa.all4$Region, summary)## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01639 0.00000 1.00000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.000000 0.000000 0.006494 0.000000 1.000000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 1.964 3.000 7.000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2447 0.0000 6.0000# no SO in LI,
tapply(sa.all4$Shrub_Oak, sa.all4$Treat_Type, summary) ## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.7094 0.0000 7.0000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.4902 0.0000 6.0000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01786 0.00000 1.00000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0# no SO in MowRX or SpringRx
rm(sa.all4)looking at SO w/o spring rx or mowrx and w/ standardized variables
sa.all_noSMRX <- sa.all3 %>%
filter(Treat_Type != "SpringRx") %>%
filter(Treat_Type != "MowRx")There is just not enough information in exposed mineral soil (so many 0s) for it to be useful. It messes up models, making them rank deficient. So just dropping that variable
so.sa1s <- glmmTMB(Shrub_Oak ~ Treat_Type + piri.ba_s + avgld_s + vegt_s + piri_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
# won't converge
#had to remove piri ba for model to converge
so.sa2s <- glmmTMB(Shrub_Oak ~ Treat_Type + avgld_s + vegt_s + piri_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC (so.sa2s) #333## [1] 333.03#avgld
so.sa3s <- glmmTMB(Shrub_Oak ~ Treat_Type + vegt_s + piri_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC (so.sa3s) #331## [1] 331.278lrtest(so.sa2s, so.sa3s) #0.6## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + avgld_s + vegt_s + piri_s + other_s +
## ba.ha_s + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + piri_s + other_s + ba.ha_s +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 16 -150.51
## 2 15 -150.64 -1 0.248 0.6185#piri
so.sa4s <- glmmTMB(Shrub_Oak ~ Treat_Type + vegt_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC (so.sa4s) #329## [1] 329.3011lrtest(so.sa3s, so.sa4s) #0.9## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + vegt_s + piri_s + other_s + ba.ha_s +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + other_s + ba.ha_s + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 15 -150.64
## 2 14 -150.65 -1 0.0231 0.8793#other
so.sa5s <- glmmTMB(Shrub_Oak ~ Treat_Type + vegt_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC (so.sa5s) #327.3## [1] 327.3061lrtest(so.sa4s, so.sa5s) #0.9 - drop other## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + vegt_s + other_s + ba.ha_s + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + ba.ha_s + offset(l.TFT) + (1 |
## Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 14 -150.65
## 2 13 -150.65 -1 0.005 0.9436#test BA
so.sa6s <- glmmTMB(Shrub_Oak ~ Treat_Type + vegt_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC (so.sa6s) #325.3## [1] 325.3149lrtest(so.sa6s, so.sa5s) #0.9## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + vegt_s + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + ba.ha_s + offset(l.TFT) + (1 |
## Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 12 -150.66
## 2 13 -150.65 1 0.0088 0.9253# test veg
so.sa7s <- glmmTMB(Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC(so.sa7s) #323.8## [1] 323.7947lrtest(so.sa7s, so.sa6s) #0.5## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -150.90
## 2 12 -150.66 1 0.4798 0.4885# add piri ba back in
so.sa8s <- glmmTMB(Shrub_Oak ~ Treat_Type + piri.ba_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC(so.sa8s) #325## [1] 325.0044lrtest(so.sa7s, so.sa8s) #0.4## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + piri.ba_s + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -150.9
## 2 12 -150.5 1 0.7902 0.374rm(so.sa1s,
so.sa2s,
so.sa3s,
so.sa4s,
so.sa5s,
so.sa6s,
so.sa8s)Best model (no springrx or mowrx)
so.sa7s <- glmmTMB(Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
summary(so.sa7s) #323.8## Family: nbinom1 ( log )
## Formula: Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## Zero inflation: ~Region
## Data: sa.all_noSMRX
##
## AIC BIC logLik deviance df.resid
## 323.8 363.6 -150.9 301.8 264
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 0.000000007557 0.00008693
## Site (Intercept) 0.109978958497 0.33163076
## Number of obs: 275, groups: Plot_No:Site, 275; Site, 26
##
## Dispersion parameter for nbinom1 family (): 0.773
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6952 0.2583 -10.435 <0.0000000000000002 ***
## Treat_TypeFallRx 0.8099 0.3689 2.195 0.0281 *
## Treat_TypeHarvest 1.5834 1.2203 1.298 0.1944
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.9034 0.7926 2.402 0.0163 *
## RegionLI 30.0378 1411573.3749 0.000 1.0000
## RegionMA 2.0029 1.3083 1.531 0.1258
## RegionME -31.2661 873734.0264 0.000 1.0000
## RegionNH -1.0435 0.8939 -1.167 0.2431
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1so.sa7s_sr <- simulateResiduals(so.sa7s, n = 1000, plot = T) #passes
testResiduals(so.sa7s_sr) #passes
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.032251, p-value = 0.9372
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9179, p-value = 0.902
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 275, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.005545455
## sample estimates:
## outlier frequency (expected: 0.00112727272727273 )
## 0## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.032251, p-value = 0.9372
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9179, p-value = 0.902
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 275, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.005545455
## sample estimates:
## outlier frequency (expected: 0.00112727272727273 )
## 0Pairwise
emmeans(so.sa7s, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.722 0.187 Inf 0.435 1.20
## FallRx 1.623 0.590 Inf 0.796 3.31
## Harvest 3.518 4.244 Inf 0.331 37.42
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 0.445 0.164 Inf 1 -2.195 0.0719
## Control / Harvest 0.205 0.250 Inf 1 -1.298 0.3965
## FallRx / Harvest 0.461 0.570 Inf 1 -0.626 0.8056
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log scalenone significantly different
Control vs FallRx close at 0.07
Now graph
Need to figure out better use of emmeans (potentially)
so.sa1.plot <- emmeans(so.sa7s, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
plot <- plot(so.sa1.plot, horizontal=FALSE, comparisons=TRUE, color = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15")) + theme_bw()
plot+
theme_classic()+
theme(panel.background = element_blank()) +
theme(axis.text = element_text(size = 14), axis.title = element_text(size = 20), legend.position = "right", legend.text = element_text(size = 14), legend.title = element_text(size = 18)) +
labs(x = "Shrub Oak Stems \n (corrected for time from treatment)",
z = "Treatment Type",
color = "Treatment Type")
More in harvest - does fire kill large scrub oak saplings? They are entirely missing from springrx and mowrx treatments
this is the same analysis, but with all treatment categories
There are no shrub oak sapling in springrx and mowrx treatments, so the above model is better - but just to have it.
so.all1 <- glmmTMB(Shrub_Oak ~ Treat_Type + piri_s + other_s + min_s + avgld_s+ piri.ba_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; function evaluation limit reached without convergence (9). See
## vignette('troubleshooting'), help('diagnose')#won't converge
#drop min & piri ba
so.all2 <- glmmTMB(Shrub_Oak ~ Treat_Type + piri_s + other_s + avgld_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.all2) #335## [1] 335.3884so.all3 <- glmmTMB(Shrub_Oak ~ Treat_Type + other_s + avgld_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.all3) #333.5## [1] 333.478so.all4 <- glmmTMB(Shrub_Oak ~ Treat_Type + avgld_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.all4) #331.6## [1] 331.5949so.all5 <- glmmTMB(Shrub_Oak ~ Treat_Type + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.all5) #329.7## [1] 329.725so.all6 <- glmmTMB(Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.all6) #327.8## [1] 327.7947so.all7 <- glmmTMB(Shrub_Oak ~ Treat_Type + piri.ba_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.all7) #329## [1] 329.0044rm(so.all1,
so.all2,
so.all3,
so.all4,
so.all5,
so.all7)
so.all6_sr <- simulateResiduals(so.all6, n = 1000, plot = T) #passes
testResiduals(so.all6_sr) #passes
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.026629, p-value = 0.8718
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.88884, p-value = 0.844
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000522088353413655 )
## 0## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.026629, p-value = 0.8718
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.88884, p-value = 0.844
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000522088353413655 )
## 0emmeans(so.all6, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.410 0.1059479 Inf 0.247 1
## FallRx 0.922 0.3350685 Inf 0.452 2
## Harvest 1.998 2.4106105 Inf 0.188 21
## MowRx 0.000 0.0000001 Inf 0.000 Inf
## SpringRx 0.000 0.0000044 Inf 0.000 Inf
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio
## Control / FallRx 0 0 Inf 1 -2.196
## Control / Harvest 0 0 Inf 1 -1.298
## Control / MowRx 10634534517938 16507512852342061056 Inf 1 0.000
## Control / SpringRx 10155053986 1101249307522083 Inf 1 0.000
## FallRx / Harvest 0 1 Inf 1 -0.626
## FallRx / MowRx 23902818311903 37103277047566393344 Inf 1 0.000
## FallRx / SpringRx 22825109079 2475234065907171 Inf 1 0.000
## Harvest / MowRx 51809117272110 80420978256787030016 Inf 1 0.000
## Harvest / SpringRx 49473193394 5365044838234219 Inf 1 0.000
## MowRx / SpringRx 0 1486 Inf 1 0.000
## p.value
## 0.1813
## 0.6927
## 1.0000
## 1.0000
## 0.9709
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
##
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scale